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SUMMARY 


Inlet flow fields for airbreathing missiles are calculated 
by the adaptation of a two-dimensional computational method 
developed for the flow around airfoils. A supersonic free 
stream is assumed to allow the forebody calculation to be 
uncoupled from the inlet calculation. The inlet calculation 
employs an implicit, time-marching finite-difference procedure 
to solve the thin-layer Navier-Stokes equations formulated in 
body-fitted coordinates. The mathematical formulation of the 
problem and the solution algorithm are given. Numerical 
stability and accuracy as well as the initial and boundary 
conditions are discussed. Instructions for program use and 
operation along with the overall program logic are also given. 


1 . 


INTRODUCTION 


This report consists of the technical explanation of the 
computer code that is employed to calculate two-dimensional 
(planar) inlet flow fields in a supersonic free stream including 
viscous effects. This code is a modified version of the program 
that was developed and used in reference 1 to calculate external 
flows over two-dimensional airfoils. In earlier reports 
(refs. 2 and 3), a detailed discussion of the inviscid version 
of this code (the Euler code) and the details of the grid-gener- 
ation routine are given. In reference 4, preliminary applica- 
tions of the method to inviscid flow in two-dimensional inlets 
are described. 

In the following sections of this report, the version of 
the code including viscous effects is explained in detail; the 
mathematical formulation of the program is given and the compu- 
tational algorithm as well as program logic are explained. A 
listing of the program and the output for a test case are 
contained in reference 5. Applications of this code and the 
Euler code to an inlet representative of current design 
practice and an assessment of the utility of the codes for such 
use (including convergence and stability behavior) are 
described in reference 6. 

2. FORMULATION OF THE PROBLEM 

In this section the governing equations (Navier-Stokes 
equations with the thin-layer approximation) for flow fields in 
two-dimensional inlets are given. Turbulent closure approxi- 
mations are described, and boundary and initial conditions as 
well as the numerical solution procedure and its stability are 
discussed . 
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2.1 Governing Equations 

The Navier-Stokes equations in Cartesian coordinates for 
plane two-dimensional flows can be written in conservation-law 
form in terms of non-dimensional variables for a perfect gas 
without external forces as (refs. 1,7) 
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S. = UT + VT + kP ^ (Y - 1) 

4 xy yy r ' y 


and 


p = (y - 1) [e - 0.5p(u^ + v^)] 


( 2 ) 


In these equations u and v are the velocities along the x 
and y coordinates respectively, p is the pressure, p is the 
density and e is the total energy per unit volume. The sound 
speed, a is given by a^ = Y(Y~1) “ 0.5(u^+v^)^ and Pr is the 

Prandtl number. Note that according to Stokes' hypothesis 
A = (-2/3) y. The reference quantities are arbitrarily selected 
but the Reynolds number (Re) is defined in terms of these 
quantities. In order to use a body-fitted coordinate system, 
the governing equations are rewritten subject to the general 
transformation . 


C = 5(x,y,t) 

n = n (x,y , t) (3) 

T = t 


where ^ is the coordinate along the body and n is the coordinate 
perpendicular to the body. In addition to this transformation, 
if the governing equations are simplified with the "thin-layer" 
approximation (ref. 8), the following equation is obtained 

9 q + 9^E + 9 F = r“^9 S (4) 

^ n e Ti 

where 

q = q/J , E = + 5 F)/J 

F = (n^q + n^E + riyF)/j 
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Here J is the transformation Jacobian given as 


J=^n - = l/(x^y - X y^) 

^x 'y ^y X ' T] y ^ 
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whereas E etc., are the metrics formed from the derivatives 

t X 

of x^ , Xp, etc., using the relations 
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Using these velocities, E and F can be written as 
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The inlet flow-field solution is obtained as the steady- 
state solution of the time-marching method employed to solve 
equation ( 4 ) . 

It should be noted that the thin-layer approximation 
neglects diffusion along the coordinate parallel to the body 
surface; all three momentum equations are retained and no 
assumption is made about the pressure. This feature enables 
the use of a single set of equations in attached and separated 
flow situations. 


2.2 The Turbulence Model 

Practical operational regimes for air-breathing missiles 
necessitate the use of a closure model to account for effects 
of transition and turbulence. In this code the thin-layer 
eddy-viscosity model of reference 8 is used. This model is a 
modified version of the two-layer eddy-viscosity model which 
was first proposed in reference 9. The modifications, which 
remove the necessity of locating the boundary-layer edge, 
consist of using the vorticity distribution to determine the 
characteristic length scale of the flow field. The Prandtl- 

van Driest formulation is used to damp the eddy-viscosity in the 
laminar sublayer. 

The basic principle of the two-layer eddy-viscosity model 
is to simulate the effects of turbulence in terms of an eddy- 




viscosity coefficient, Hence the molecular viscosity 

coefficient, y , that appears in the Navier-Stokes equations is 
replaced by y + and y^ is calculated separately for the 

inner and outer regions of the boundary layer, i.e. 
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where y is the normal distance from the wall and y 

^ -^crossover 

is the smallest value of y at which values from the inner and 
outer formulas are equal. 

The Prandtl-Van Dreist formulation is used in the inner 
region to damp y^, 

(y^) = |o3 I (10) 

inner 


where i (the turbulence length scale) is defined as 


£ = ky[l - exp (-y’*’/A'^) ] (11) 

Here 1 w | is the magnitude of the vorticity vector and for 
two-dimensional flows is given as 
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and 
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The outer layer eddy-viscosity coefficient is calculated 


from 


(y 
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( 14 ) 


where K is the Clauser constant, C^p is an additional constant, 
and 
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The quantities and are determined from the function 

^ MAX MAX 


F(y) 


y CO [1 


exp (-y'^/A"'’) ] 


( 16 ) 


The quantity is the maximum value of F (y) that occurs in 
a profile and Y^yj^^ is the value of y at which it occurs. The 
function F (y) is the Klebanoff intermittency factor given 
by 
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The quantity is the difference between maximum and minimum 

L) Xr 

total velocity in the profile (i.e., at a fixed x station) 
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The second term in is taken to be zero (except in wakes) . 

The outer formulation (eqs. (14) and (15)) can be used in 
both attached and separated boundary layers. 

In this model the effect of transition to turbulence is 
simulated by setting equal to zero everywhere in a profile 
for which the maximum tentatively computed value of y^ from 
the foregoing relations is less than a specified value, that is. 


y^ = 0 if 



max in 
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The constants appearing in the foregoing relations have 
been determined by requiring agreement with the formulation 
of reference 8 for constant pressure boundary layers at tran- 
sonic speeds. The values thus determined are 
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2 . 3 Boundary and Initial Conditions 


In the problem under consideration there are four types 
of boundary conditions (fig. 1). These are solid boundaries, 
inflow boundaries, outer boundaries and outflow boundaries. 

Along the cowl surface and the ramp surface the no-slip 
viscous boundary condition is satisified, i.e. U = 0 and V = 0. 
The pressure on the body surface is found from the normal 
momentum equation. This gives 


(n C +n C ) + P (n^+n^) = o 

^ 'x X y y n x y 


( 20 ) 


Throughout the integration procedure values of p at the 

n+1 

advanced time step (p ) are assumed equal to values of p at 
the current time step (p^) on the flow-field boundaries; the 
latter are obtained by linear extrapolation from the interior. 
This simplifies the procedure for implementing the boundary 
conditions but results in first-order accuracy in time on the 
boundaries. 

Free stream values are specified at the inflow and outer 
boundaries. These may be non-uniform values as obtained from 
a solution for a forebody as described in reference 5. In the 
case of a solution including the effects of a forebody, the 
outer boundary is extended such that it contains the shock from 
the forebody. For computational economy in the use of mesh 
points in the work done to date, the inflow boundary has been 
located a short distance downstream of the leading edge of the 
ramp. This requires special treatment for the mesh points on 
the inflow boundary which are downstream of the ramp-leading- 
edge shock. This treatment is discussed below under the 
paragraph dealing with initial conditions. 

The outflow boundary conditions on the portion of the 
outflow boundary external to the inlet are calculated by zeroth- 
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order extrapolation from the interior. On that portion inside 
the inlet duct, the boundary conditions are specified according 
to whether the condition being calculated is supercritical or 
subcritical. In the case of supercritical operation, the flow 
field variables are again calculated by zeroth-order extrapola- 
tion from the interior. For subcritical operation (subsonic 
outflow) a set of boundary conditions are incorporated which 
are obtained from the steady-state forms of the governing 
equations assuming uniform parallel outflow. To allow the use 
of this assumption, a constant-area section is added to the 
downstream portion of the inlet duct such that v = 0. Back 
pressure is prescribed consistent with subsonic outflow (constant 
along y-direction) , and u and p are found by zeroth-order 
extrapolation from the interior. The energy, e, is calculated 
from the values of u, p and p. We have found em.pirically that 
in order to conserve total enthalpy for a subcritical calcu- 
lation, the high pressure boundary condition has to be 
introduced over some 300 time steps. A branching is also 
provided to test whether the flow is subsonic or supersonic 
at the outflow boundary (the gradual increase in pressure, 
of course, will make the flovz subsonic gradually). If at a 
particular time step the flow is supersonic, all the variables 
are found by zeroth-order extrapolation, otherwise the above 
subsonic conditions are imposed. 

The initial conditions are specified by using either the 
impulsive initial condition, i.e., free-stream conditions are 
imposed throughout the flow field, or the final solution of a 
previously calculated flow field is used. When starting 
from free-stream conditions, inviscid values calculated from 
2-D shock theory are used at the 3-5 points on the inflow 
boundary that are downstream of the ramp shock. No allowance 
is made for the ramp boundary layer at the inflow boundary. 
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2 . 4 Solution Procedure and the Numerical Schemes 


In the solution procedure employed to integrate equation (4) 
a temporal linearization process is employed. In order to 
employ this on the vector terms of equation (4) one needs to 

/\ ^ /S /N ^ 

evaluate the Jacobian Matrices A e 3E/3q and B e 3F/3q. 

/N 

The flux vectors E and F are both linear combinations of 
q, E, and F (ref. 1) . The viscous term is linearized by using 
the principle of homogeneity (ref. 1) such that 


gn+l ^ gn ^ j-1 -n+1 


( 21 ) 


where M is a 4x4 matrix whose elements are functions of q. 

In the solution procedure , after approximate factorization 
the linearized equations are cast into delta-form, (AF) algorithm 
(refs. 1, 10). This algorithm is non-iterative, and requires 
the inversion of two-block-tridiagonal (4x4) coefficient 
matrices at each time step in the integration procedure at 
only the interior points. References 1, 4, 10, 11 and 12 contain 
detailed descriptions and various applications of the delta- 
form AF algorithm, hence here we outline it only briefly. The 
delta-form AF algorithm can be used either with trapezoidal 
or Euler temporal implicit differencing and reads 

(I+h6 A^) (I+h6 B^+h6 M^) (q’^'^^-q^) 

K 0 n 

= -At(6.E^+6 F^-r“^ 5 S^) - aAt(6-^^^+6 ) (22) 

? n e n KKKK nrinn ^ 

Here 6^ and 6^ are second-order central-difference operators. 

At is the integration step size, and h = At or At/2 for first- 
order or second order two-level time differencing, respectively. 
The fourth-order smoothing terms are added to the RHS of the 
difference equation both to overcome non-linear instability 
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and to damp short wave lengths (refs. 1, 10 and 11). It is 

generally desirable to add the smoothing term in terms of 

both explicit and implicit portions. The explicit part 

suppresses non-linear instabilities whereas the implicit part 

enables the use of a ~ 0 (At) for very large At where a is the 

artificial viscosity coefficient. In the current work, 

these terms have been included in a pseudo-implicit manner by 

1 

inverting a product of scalar pentadiagonals 


-aAt ( 6 
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It should be noted that for compatibility with the 
boundary conditions, a second-order smoothing is used at points 
adjacent to the boundary, whereas no smoothing is used at the 
boundary points themselves. 


2.5 Numerical Stability 

In explicit methods, the step size in the marching 
direction is bounded by the Courant-Fr iedr ics-Lewy (CFL) 
condition. According to this, for block-tridiagonal matrix 
inversion along the AF scheme requires 


At 

A^ 


A max 


< 1 


(24a) 


where (cr^) is the maximum spectral radius of the local eigen- 

riia.X ^ 

/N 

values of A and the expression given by equation (24a) is 
referred to as the Courant number. Similarly, for block- 
tridiagonal matrix inversion along n, one can write the Courant 


13 



■IIHIIIIHIIIIII I 


il I 


I II I II I 


l_ 


II III II 


number as 


At 

An 


B max 


< 1 


(24b) 


where (cr^) is the maximum spectral radius of the local eigen- 
B max 

values of B. A practical constraint in a two-dimensional 
problem can be written as 


At max of ^ 1 


On the basis of linear stability analysis, in implicit schemes 
this condition can be relaxed so that the right-hand-side 
constant can assume values much larger than one. In the present 
work, however, we have found that this was possible only for 
the case of supercritical operation with inviscid flow. In 
this case typical values were around 5 to 10. In the cases of 
inviscid flow with subcritical outflow boundary conditions 
and for viscous flow, for most of the calculation the Courant 
number had to be kept smaller than one to achieve numerical 
stability . 

It should also be noted, however, that in the AF algorithm, 
the modulus of the amplification matrix ( |X|, ref. 1) is close to 
one regardless of the magnitude of At and therefore the scheme 
is only neutrally stable {for a stable scheme |a| should 
decrease with increasing At). Therefore, steady-state 
convergence cannot be accelerated by the use of very large 
time-steps, even where their use introduces no instability. 
Further information on the stability and convergence of this 
code can be found in reference 6. 

With the inclusion of viscous terms, additional restrictions 
on the stability of the numerical scheme are required. Usually 
referred to as the "diffusive stability" restrictions, these 
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involve the prescription of a mesh system that ensures Re^ - 2, 
where Re^ is the mesh Reynolds number. In high Reynolds number 
flows with Rql 1 (where L is a characteristic body dimension) , 
this criterion will usually be met if there is fine resolution in 
the cross-stream, y-direction, e.g. Ay - 6/10 (ref. 13). A coarse 
mesh in the x-direction, Ax - L/10 is generally accepted as ade- 
quate. In the present viscous flow calculation (Re = 3.28 x 10^/ 
meter) grid spacing in the n~cii^ection varied exponentially with 
distance from the wall. The minimum spacing close to the cowl and 
ramp boundaries was around 0.3048 x 10~^ meters (<< 6/10) in accord 
with the requirements of the turbulence model; so the requirements 
for diffusive stability were fulfilled. 

3. OVERALL PROGRAM LOGIC 


In this section the overall logic of the computer program 
used to solve equation (22) is summarized. The details are 
given in reference 1. The computational procedure starts by 
forming the right-hand side terms (RHS) in equation (22) , first 
the smoothing operator then the steady part. The computational 
algorithm there forms Aq by block tridiagonal matrix inversion 
in ^ 


(I + h6^A^) Aq* = RHS (26) 

^ . "'n+l 

Aq* is stored temporarily in q . The procedure continues 

with block tridiagonal matrix inversion in ri 


(I + h6^B^) Aq = Aq- 


^ ri"i" 1 

Finally q is obtained from 


^n+1 "n ."n 
q = q + Aq 


(27) 


(28) 
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Note that for the two-level scheme used here, the algorithm 
requires tvjo levels of data storage. 

In the computer program used in this work the execution 
of the above outlined algorithm is performed by PROGRAM MAIN. 
It starts the computation by calling SUB.INITIA. This 
subroutine defines the physical and mathematical constants 
that are to be used in the computational procedure. The 
initial field is also specified in this subprogram. An option 
is provided for starting the calculation from a previously 
calculated (converged) field. This subroutine also calls 
SUB, GRID by which either a simple-stretched grid for a flat 
plate is computed or a previously calculated grid is read from 
a disc. Then SUB . XYMETS is called to evaluate the metric 
coefficients for the grid system. 

The main program then calls SUB. EIGEN, v/hich calculates 

/N /\ 

the eigenvalues of the coefficient matrices A and B; SUB. EIGEN 
also finds the spectral radius and finally computes the time- 
step from the prescribed Courant number. The optional SUB. MAP 
maps the P-matrix; this can either be the Jacobian of trans- 
formation or the calculated pressure field. The normalized 
variables and P/P^ of the initial field are then printed by 
SUB. OUTPUT. Free-stream quantities are given or computed by 
SUB.HARVIO. 

After the initial field is prescribed and printed, the 
main program starts to execute the integration loop. This is 
done by calling SUB. STEP. This subroutine in turn calls 
SUP.BC to calculate data at the flow field boundaries. 

SUB. STEP then calls SUB.RHS which computes the steady part of 
the difference algorithm according to whether the flow is 
inviscid or viscous. If the variable INVIS is set equal to or 
less than zero, the program does not account for viscous 
effects. If INVIS > 0, SUB.VISRHS is called. The specified 


16 



variable LAMIN determines if laminar or turbulent viscous 

effects are to be considered. If LAMIN £ 0, SUB.MUKIN is called 

to calculate the molecular viscosity coefficient, p, from 

Sutherland's formula. If LAMIN > 0, SUB.MUTUR is called to 

calculate the eddy viscosity coefficient and SUB.MUKIN is 

called to calculate y. Finally, the effective viscosity 

y = y^. + y is formed. SUB. STEP then calls SUB. SMOOTH which 
eft t 

forms the smoothing operator. The temporal differences are 
evaluated by SUB. DIFFER. SUB. STEP then calls SUB.FILTRX which 
performs the block tridiagonal matrix inversion in the 
^-direction and SUB.FILTRY which performs the block tridiagonal 
matrix inversion in the y-direction . Final output is furnished 
by subroutine SUB.HARVIO at the end of the iteration loop. 

The original program contains options for upwind differencing 
in conservative form (ref. 1) to be used in transonic flow 
calculations. These subroutines have not been used in the 
present work. 


4 . PROGRAM USE AND OPERATION 

In this section we give a detailed description of the 
input information required to run a typical case. As described 
m Section 3 , the initial field as well as values of various 
constants are evaluated in SUB.INITIA. This subroutine also 
reads the input data except for the last card, which specifies 
geometry (planar or axisymmetr ic ) and type of outflow boundary. 
The last card is read by SUB.BC. We now list the input data 
cards in order: 



Card No. 


Format 


Variables 


1 


815 


NMAX: 


jriAX: 


KMAX: 

NP: 


815 


METH: 


I RE AD: 
INVIS : 


2 


715 


IREGO: 


3 


ISTORE: 

JTAILl : 
JTAIL2 : 


315 


IPLOT: 
I UPWIND: 


lOSCIL: 

LAMIN: 


4 


8F10. 0 


CNBR: 

DX: 


Maximum number of time steps 
Maximum number of points along ^ 
Maximum number of points along n 
Number of time steps for calling 
SUB . MAP . 

Flag for upwind differencing. If 
METH > 0, program skips upwind 
differencing . 

Option for grid system. If IREAD 
> 0, reads grid from disc. 

If INVIS > 0, program accounts for 
viscous effects; if INVIS = 0, the 
calculation is inviscid. 

If IREGO > 0, reads initial field 
from disc, otherwise starts from 
free-stream conditions. 

If ISTORE > 0, solution data 
stored on disc. 

First ^-interior point, set to 1. 
Last ^-interior point, set to 
JMAX-1 . 

Index for calling plot routine. 

If lUPWIND > 0, skips upwind 
differencing . 

For stationary airfoil set 
lOSCIL = 0 

If LAMIN > 0, calculates turbulent 
viscosity . 

Courant number; set to about 10. 
x-increment used in SUB. GRID. For 
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Format 


Variables 


Card No. 


grid read off the disc set DX = 0 . 
DY : y-increment, as above. 

FSMACH: Free-stream Mach number 

SMU : Pseudo viscosity coefficient. Set 
to about 10 times the time step. 
8F10.0 EPS: Used in SUB. GRID. Here set 

EPS = 0. 

RE: Reynolds number 
ALPHA; Angle of attack 

4F10.0 I XOSCIL: | Variables required for oscil- 
VARA: I lating airfoil. For the inlet 
VARB: ( flow set to zero. 

VARC: ) 


6 


212 


LFAC : If LFAC = 1 outflow is subsonic, 

otherwise it is supersonic 
JAXI : If JAXI = 0 flow is plane 2-D, 

If JAXI = 1 flow is axi symmetric 


(Note that axisymmetric option has been added only for the 
inviscid case. When calculating viscous flows, set JAXI = 0) . 

In starting from free-stream conditions, the number of 
points on the inflow boundary to be placed downstream of the 
ramp shock are prescribed in SUB.INITIA. Values for the 
variables are also prescribed in this subroutine using 2-D 
oblique shock theory. 

The typical output of the program involves the following 
(described in more detail in reference 5) . 

(1) Tabulation of the input parameters. 

(2) Printout of the Jacobian matrix as calculated from 
the generated grid. 
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(3) A map of the transformation Jacobian. 

(4) Printout of the flow field variables on the cowl 
surface. 

(5) Free -stream quantities given by SUB.HARVIO. 

(6) Iteration index, physical time and time step. 

(7) Maximum ERRORS occurring within the field. 

(8) Values of LFAC (if LFAC > 1 outflow is subsonic) and 
JAXI. 

(9) Mass flow ratio evaluated at various ^-stations. 

(10) Residual quantity; must converge to a small number such 
as O(10"^) - O(10”^) . 

(11) Line-printer plot of pressure in the computational 
plane . 

(12) Line-printer plot of Cp on the ramp surface. 

(13) Line-printer plot of Cp on the cowl surface. 

(14) At time step equal to NMAX total output as furnished 
by HARVIO. 

Note that the program that has been employed to generate 
the body-fitted coordinate system has been documented else- 
where (refs. 2,3). 

5. PROGRAM NUMERICAL ACCURACY AND LIMITATIONS 

The local accuracy of the finite difference scheme that is 
employed in this program is second-order both in time and space. 
It should be noted that the method is a conservative difference 
formulation, therefore its global numerical accuracy is also 
given by the local accuracy and is second-order in time and 
space. There are, however, other points that contribute to 
the accuracy of the numerical solution: 

(a) The grid system generated has a very important effect 
on the overall solution. With a course grid, the shocks 
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(discontinuities) are not adequately resolved, and in certain 
cases the resulting field is totally erroneous. Near the solid 
walls the grid system must be properly clustered. Grids also 
should vary smoothly in the flow field. 

(b) The overall accuracy also depends on the accuracy at 

the boundaries. In the present program at the flow-field 

. . ^n+1 ^n ^ 

boundaries xt xs assumed that q equals q . Therefore on the 

boundaries the formal accuracy is first-order in time. Second 

order accuracy in space, however, is retained all through the 

flow field. 

The major formal limitations of the program are (a) the 
flow field analyzed is two-dimensional, and (b) owing to the 
neutral stability of the method, convergence to steady state 
cannot be accelerated by using large time-steps. The accuracy 
of the results obtained in a specific application depends on 
the grid, as explained above, and must be evaluated ultimately 
by comparison with experimental data. Also the adequacy of the 
turbulence model used has not been addressed and requires 
separate study. 
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